Pour cette deuxième séance, nous allons apprendre à utiliser un planifateur d’itinéraire sur R afin de calculer des isochrones autour de chaque supermarché et ainsi pouvoir identifier les zones de chalandises.
Une première information simple à obtenir mais très informative peut être obtenue en traçant des courbes d’isodistances autour de chaque supermarchés. On peut ainsi indentifier les zones non-deserviés.
télécharger le package en .zip en cliquant sur ce lien
https://cran.r-project.org/src/contrib/Archive/opentripplanner/opentripplanner_0.4.0.tar.gz
Et installer Java 8 sur votre ordinateur:
https://www.java.com/fr/download/
library(dplyr)
library(ggplot2)
library(sf)
library(readr)
library(readxl)
library(tmap)
library(tmaptools)
library(here)
otp_available <- requireNamespace("opentripplanner", quietly = TRUE)
if (!otp_available) {
otp_pkg <- here("Cours", "opentripplanner_0.4.0.tar.gz")
if (file.exists(otp_pkg)) {
try(install.packages(otp_pkg, repos = NULL, type = "source"), silent = TRUE)
otp_available <- requireNamespace("opentripplanner", quietly = TRUE)
}
}
## Warning in install.packages(otp_pkg, repos = NULL, type = "source"): 'lib =
## "/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library"' is not
## writable
if (otp_available) {
library(opentripplanner)
} else {
message("Package 'opentripplanner' absent: les sections OTP/isochrones seront ignorees.")
}
# Chargement des données locales (mêmes sources que test.Rmd)
iris_shp <- st_read(
dsn = here("Cours"),
layer = "georef-france-iris-millesime"
)
## Reading layer `georef-france-iris-millesime' from data source
## `/Users/loic/Documents/UGA/donnee_entreprise/ProjetECO/Cours'
## using driver `ESRI Shapefile'
## Simple feature collection with 728 features and 30 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4.74204 ymin: 44.69597 xmax: 6.359009 ymax: 45.88339
## Geodetic CRS: WGS 84
sirene_shp <- st_read(
dsn = here("Cours"),
layer = "dataseed-sirene-2"
)
## Reading layer `dataseed-sirene-2' from data source
## `/Users/loic/Documents/UGA/donnee_entreprise/ProjetECO/Cours'
## using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 342 features and 106 fields (with 5 geometries empty)
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 4.767436 ymin: 44.90643 xmax: 6.121839 ymax: 45.81497
## Geodetic CRS: WGS 84
iris_data <- read_excel(here("Cours", "iris_isere.xlsx"))
# Zone d'étude: Grenoble + communes nord
sirene_coms <- c(
"GRENOBLE", "MEYLAN", "LA TRONCHE",
"CORENC", "SAINT-MARTIN-LE-VINOUX", "SAINT-EGREVE"
)
iris_coms <- c(
"Grenoble", "Meylan", "La Tronche",
"Corenc", "Saint-Martin-le-Vinoux", "Saint-Égrève"
)
sirene_grenoble <- sirene_shp %>%
filter(libellecomm %in% sirene_coms & etatadminis == "Actif")
iris_grenoble <- iris_shp %>%
filter(com_name %in% iris_coms)
# Préparation des variables IRIS (même logique que test.Rmd)
colnames(iris_grenoble)[20] <- "CODE_IRIS"
iris_data$nb_residence_princ <- iris_data$Res_princ_30_moins_40_m2 +
iris_data$Res_princ_40_moins_60_m2 +
iris_data$Res_princ_60_moins_80_m2 +
iris_data$Res_princ_80_moins_100_m2 +
iris_data$Res_princ_100_moins_120_m2
vars_to_keep <- c(
"CODE_IRIS",
"1er_quartile_euro", "3e_quartile_euro",
"Actifs_15-64_ans2", "Pop_15-64_ans",
"Delta_Pop_20-64_ans", "nb_residence_princ"
)
iris_grenoble <- merge(
iris_grenoble,
iris_data[, vars_to_keep],
by = "CODE_IRIS",
all.x = TRUE
)
iris_grenoble <- iris_grenoble %>%
rename(
premier_quartile_revenu = `1er_quartile_euro`,
troisieme_quartile_revenu = `3e_quartile_euro`,
Actif_total = `Actifs_15-64_ans2`,
Pop_15_64 = `Pop_15-64_ans`,
Delta_pop_20_64ans = `Delta_Pop_20-64_ans`,
Nb_residence_principale = nb_residence_princ
)
iris_grenoble <- iris_grenoble %>%
mutate(
across(
c(
premier_quartile_revenu,
troisieme_quartile_revenu,
Actif_total,
Pop_15_64,
Delta_pop_20_64ans,
Nb_residence_principale
),
~ na_if(., 0)
)
)
st_crs(sirene_grenoble)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
# Buffer de 500 m autour des supermarchés
sirene_grenoble_l93 <- st_transform(sirene_grenoble, 2154)
buffer <- st_buffer(sirene_grenoble_l93, 500)
buffer <- st_transform(buffer, st_crs(iris_grenoble))
tmap_mode("view")
tm_shape(iris_grenoble) +
tm_text(text = "iris_name", size = 0.8) +
tm_polygons(
"Pop_15_64",
alpha = 0.5,
style = "pretty",
id = "iris_name",
title = "Population (15-64 ans)"
) +
tm_shape(buffer) +
tm_borders("chartreuse3", lwd = 3) +
tm_shape(sirene_grenoble) +
tm_dots(
"enseigne1et",
legend.show = TRUE,
title = "Enseigne",
id = "enseigne1et",
popup.vars = "trancheeffe.2",
alpha = 0.3
) +
tm_layout(legend.outside = TRUE)
tm_shape(iris_grenoble) +
tm_text(text = "iris_name", size = 0.8) +
tm_polygons(
"Pop_15_64",
alpha = 0.5,
style = "pretty",
id = "iris_name",
title = "Population (15-64 ans)"
) +
tm_shape(buffer) +
tm_polygons("enseigne1et", legend.show = FALSE, id = "enseigne1et", alpha = 0.5) +
tm_shape(sirene_grenoble) +
tm_dots(
"enseigne1et",
legend.show = TRUE,
title = "Enseigne",
id = "enseigne1et",
popup.vars = "trancheeffe.2",
alpha = 0.8
) +
tm_layout(legend.outside = TRUE)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "pretty"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: use `fill.legend = tm_legend_hide()` instead of
## `legend.show = FALSE`.
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
## [tm_dots()] Arguments `legend.show` and `title` unknown.
Un dossier & un sous-dossier specifiques à créer dans votre R working directory > OTP > graphs > default
Des données sur le réseau de transport GTFS (General Transit Feed Specification) à déposer dans OTP > graphs > default
Un fond de carte openstreetmap à déposer dans OTP > graphs > default
gtfs_img_candidates <- c(
here("Cours", "gtfs.png"),
here("Cours", "Capture.PNG"),
here("Cours", "Capture.png")
)
gtfs_img <- gtfs_img_candidates[file.exists(gtfs_img_candidates)][1]
if (!is.na(gtfs_img)) {
knitr::include_graphics(gtfs_img)
} else {
cat("Ajoutez une capture GTFS dans Cours/gtfs.png pour l'affichage.")
}
Vous pouvez par exemple trouver les données GTFS sur le site M mobilité du TAG: data.mobilites-m.fr
Une fois téléchargé le fichier SEM-GTFS.zip renommez le en gtfs.zip et placez le dans OTP > graphs > default
Il y a aussi beaucoup de sources, par exemple :
https://download.geofabrik.de/europe/france/rhone-alpes.html
Placez le fichier dans OTP > graphs > default
path_data <- NULL
for (candidate in c(here("Cours", "OTP"), here("Cours", "data"))) {
if (
dir.exists(file.path(candidate, "graphs", "default")) ||
file.exists(file.path(candidate, "Graph.obj"))
) {
path_data <- candidate
break
}
}
if (is.null(path_data)) {
path_data <- here("Cours", "OTP")
dir.create(path_data, recursive = TRUE, showWarnings = FALSE)
}
path_otp <- otp_dl_jar()
# Vous n'avez besoin de lancer la construction du graphe qu'une fois
# log <- otp_build_graph(otp = path_otp, dir = path_data, memory = 15000)
log1 <- otp_setup(otp = path_otp, dir = path_data)
otpcon <- otp_connect()
otp_img_candidates <- c(
here("Cours", "OTP_capture.png"),
here("Cours", "OTP_capture.PNG")
)
otp_img <- otp_img_candidates[file.exists(otp_img_candidates)][1]
if (!is.na(otp_img)) {
knitr::include_graphics(otp_img)
} else {
cat("Ajoutez une capture OTP dans Cours/OTP_capture.png pour l'affichage.")
}
iso <- otp_isochrone(
otpcon,
fromPlace = sirene_grenoble,
fromID = sirene_grenoble$siret,
mode = c("WALK", "TRANSIT"),
maxWalkDistance = 2000,
date_time = as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
cutoffSec = c(5, 7, 9, 11) * 60
)
iso$minutes <- iso$time / 60
# Vérification
colnames(iso)[3] <- "siret"
iso <- merge(
iso,
st_drop_geometry(sirene_grenoble[, c("siret", "enseigne1et")]),
by = "siret",
all.y = TRUE
) %>%
unique()
# Option de correction ponctuelle (si un point ne retourne pas d'isochrone)
if (nrow(sirene_grenoble) >= 16) {
st_geometry(sirene_grenoble)[16] <- st_sfc(
st_point(c(5.728253270386639, 45.17659752397455)),
crs = st_crs(sirene_grenoble)
)
}
# On recommence
iso <- otp_isochrone(
otpcon,
fromPlace = sirene_grenoble,
fromID = sirene_grenoble$siret,
mode = c("WALK", "TRANSIT"),
maxWalkDistance = 1500,
date_time = as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
cutoffSec = c(5, 7, 9, 11) * 60
)
colnames(iso)[3] <- "siret"
iso <- merge(
iso,
st_drop_geometry(sirene_grenoble[, c("siret", "enseigne1et")]),
by = "siret",
all.y = TRUE
) %>%
unique()
iso$minutes <- iso$time / 60
tmap_mode("view")
tm_shape(iris_grenoble) +
tm_polygons(
"Pop_15_64",
alpha = 0.5,
style = "pretty",
id = "iris_name",
title = "Population (15-64 ans)"
) +
tm_shape(sirene_grenoble) +
tm_dots(
"enseigne1et",
legend.show = TRUE,
title = "Enseigne",
size = 0.1,
id = "enseigne1et",
popup.vars = "trancheeffe.2"
) +
tm_layout(legend.outside = TRUE) +
tm_shape(iso) +
tm_fill(
"minutes",
breaks = c(0, 5.01, 7.01, 9.01, 11.01),
title = "Isochrones (minutes)",
style = "fixed",
labels = c("0 to 5", "5 to 7", "7 to 9", "9 to 11"),
palette = "-BuPu",
id = "minutes",
alpha = 0.3
) +
tm_borders()
iso420 <- iso %>% filter(time == 420)
iris_grenoble$nb_supermarche_per_iris <- rowSums(
ifelse(st_intersects(x = iris_grenoble, y = iso420, sparse = FALSE), 1, 0)
)
tm_shape(iris_grenoble) +
tm_polygons(
"nb_supermarche_per_iris",
style = "pretty",
id = "iris_name",
labels = c("0", "1", "2", "3", "4", "5"),
title = "Nb de Supermarchés accessible en 7min"
) +
tm_text(text = "iris_name", size = 0.8)
data_temp <- matrix(0, nrow = nrow(iris_grenoble), ncol = nrow(iso420))
for (j in 1:nrow(iso420)) {
for (i in 1:nrow(iris_grenoble)) {
data_temp[i, j] <- as.character(st_intersects(iris_grenoble, iso420, sparse = FALSE)[i, j])
data_temp[i, j] <- ifelse(data_temp[i, j] == "TRUE", iso420$enseigne1et[j], 0)
}
}
data_temp2 <- matrix(0, nrow = nrow(data_temp), ncol = length(unique(iso420$enseigne1et)))
for (j in 1:length(unique(iso420$enseigne1et))) {
data_temp2[, j] <- apply(
data_temp,
1,
function(x) length(which(x == unique(iso420$enseigne1et)[j]))
)
}
data_temp2 <- as.data.frame(data_temp2)
colnames(data_temp2) <- unique(iso420$enseigne1et)
write.csv(data_temp2, here("Cours", "nb_pdv_iris.csv"), row.names = FALSE)
intersect_pct <- st_intersection(iris_grenoble, iso420) %>%
mutate(intersect_area = st_area(.)) %>%
dplyr::select(iris_name, intersect_area) %>%
group_by(iris_name) %>%
mutate(intersect_area = sum(intersect_area)) %>%
st_drop_geometry() %>%
unique()
iris_grenoble <- mutate(iris_grenoble, iris_area = st_area(iris_grenoble))
iris_grenoble <- merge(iris_grenoble, intersect_pct, by = "iris_name", all.x = TRUE)
iris_grenoble <- iris_grenoble %>%
mutate(intersect_area = ifelse(is.na(intersect_area), 0, intersect_area)) %>%
mutate(couverture = as.numeric(intersect_area / iris_area) * 100)
tm_shape(iris_grenoble) +
tm_polygons(
col = "couverture",
alpha = 0.5,
style = "cont",
id = "iris_name",
title = "Couverture (%)"
) +
tm_text(text = "iris_name", size = 0.8) +
tm_shape(sirene_grenoble) +
tm_dots(
"enseigne1et",
legend.show = TRUE,
title = "Enseigne",
size = 0.1,
id = "enseigne1et",
popup.vars = "trancheeffe.2"
) +
tm_layout(legend.outside = TRUE)
st_write(
iris_grenoble,
here("Cours", "iris_grenoble2.gpkg"),
delete_dsn = TRUE
)
Préparez plusieurs cartes avec des isochrones différentes et d’autres variables au niveau des IRIS (carte 3.1) afin de proposer une localisation potentielle d’un nouveau supermarché à Grenoble. Proposez aussi une autre carte (3.2) construite à partir d’une isochrone différente (>7min). Ce travail sera évalué.